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We investigate by direct numerical simulations the flow that rising bubbles cause in an originally 
quiescent fluid. We employ the Eulerian-Lagrangian method with two-way coupling and periodic 
boundary conditions. In order to be able to treat up to 288000 bubbles, the following approximations 
and simplifications had to be introduced, as done before e.g. by Climent and Magnaudet, Phys. 
Rev. Lett. 82, 4827 (1999): (i) The bubbles were treated as point-particles, thus (ii) disregarding 
the near-field interactions among them, and (iii) effective force models for the lift and the drag 
forces were used. In particular, the lift coefficient was assumed to be 1/2, independent of the 
bubble Reynolds number and the local flow field. The results suggest that large scale motions are 
generated, owing to an inverse energy cascade from the small to the large scales. However, as the 
Taylor- Reynolds number is only in the range of 1, the corresponding scaling of the energy spectrum 
with an exponent of —5/3 cannot develop over a pronounced range. In the long term, the property 
of local energy transfer, characteristic of real turbulence, is lost and the input of energy equals 
the viscous dissipation at all scales. Due to the lack of strong vortices the bubbles spread rather 
uniformly in the fiow. The mechanism for uniform spreading is as follows: Rising bubbles induce a 
velocity field behind them that acts on the following bubbles. Owing to the shear, those bubbles 
experience a lift force which make them spread to the left or right, thus preventing the formation 
of vertical bubble clusters and therefore of efficient forcing. Indeed, when the lift is artifically put 
to zero in the simulations, the fiow is forced much more efficiently and a more pronounced energy 
accumulates at large scales (due to the inverse energy cascade) is achieved. 



I. INTRODUCTION 



fore bubble injection. 



Bubbly flows are ubiquitous in nature and technology, 
but exact analytical or numerical results are extremely 
difficult to obtain, due to the multiscale nature of the 
problem and due to the moving interfaces. For a review 
on the numerical modeling we refer to ,11, [2] , for recent 
reviews on the experimental situation we refer to 0, |3| j 
and our own recent work on the subject has been summa- 
rized in [sj. An excellent overview on the experimental, 
numerical, and theoretical knowledge for various bubble 
Reynolds numbers can also be found in refs. 0, 0|- 

The motion of the small bubbles in the fluid induces 
velocity fluctuations that can be either dissipated imme- 
diately by viscosity or can be enhanced, thus generating 
motion on scales much larger than the disturbance di- 
mension. Owing to their random character, these fluctu- 
ations are refered to as "pseudo-turbulence" 0, d, Hi • In 
a flow initially at rest and only forced by rising bubbles, 
the pseudo-turbulent fluctuations are the only source of 
energy. Otherwise they can add to the already existing 
fluid velocity fluctuations, which are driven in some other 
way. Depending on the flow conditions and the bubble 
size distribution, the turbulent energy dissipation may be 
either enhanced or suppressed. 

Lance and Bataille [9] have suggested that the effect 
of bubbles on the flow depends on the ratio of kinetic 
energy due to the bubble motion and the typical kinetic 
energy (u'^) of the fluctuations in the liquid velocity be- 



(1) 



where u' is the typical flow velocity fluctuation, a the 
void fraction, vt the bubble rise velocity in still water, 
and we have taken i for the added mass coefficient. The 
ratio 6 is called the bubblance parameter For 6^1 
the kinetics of the bubbly flow are entirely dominated 
by the turbulent energy of the flow and the bubbles can 
be considered as some distortion, such as in the experi- 
ments of refs. O, [nUllJllJI^^ 



or in the numerical 



simulations of refs. [l6l IITI" [2]| . whereas for 

6^1 the flow is dominated by the bubble motion; i.e. 
we have bubblance rather than turbulence, such as in 
the experiments of ref. [1, lol. J2 ^ . l23l [2^. [25l or in the 
numerical simulations of ref. [mTl. l26l. l27l. l28l. [29j. The 
analogous situation for sedimenting particles has exper- 
imentally been analyzed by Faeth and coworkers, both 
for particles in water [s^l and in air [sj] . 

In the present work we focus on the bubblance case 
& ^ I, namely on microbubbles rising in an initially qui- 
escent flow, where formally b = oo. These conditions 
imply that pseudo-turbulence, due to the bubbles' buoy- 
ancy, is the only source of flow energy, thus bubbles drive 
the turbulence and eventually the energy dissipation. We 
will address the following questions: (i) What is the time 
evolution of the energy of bubbly driven turbulence ini- 
tially at rest? (ii) Are microbubbles able to induce in 
still fluid a flow that possesses similar features as real 
turbulence, i.e., can the inertial scaling law characteris- 
tic of homogeneous and isotropic turbulence be attained 
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in a flow forced solely by bubbles? (iii) How are bub- 
bles eventually distributed in such flow and what forces 
determine this distribution? 

To be able to address these questions, the flow should 
best be driven at least by ten-thousands of bubbles, in 
order to have sufficient statistics. While hundreds of bub- 
bles can still be treated with front-tracking techniques 
@, 0, m HE H &2|], resolving both the shape of the 
bubble and the flow around it, this clearly is no longer 
possible for ten-thousands of bubbles, and one therefore 
is forced to employ approximations. We will thus fol- 
low the complementary approach and model the bub- 
bles as point-particles, on which effective forces such as 
drag and lift act similarily as has also been done 
by Climent and Magnaudet 2u\. We thus disregard the 
near-flcld interactions among the bubbles. Also the effec- 
tive forces had to be modelled, namely by chosing drag 
and lift coefficients. Unfortunately, the lift coefficient 
is not well known In fact, it depends on the bub- 
ble Reynolds number and on the local shear and vor- 
ticity in a highly non-trivial way [H [13, [H, HI, [13, HI] . 
Moreover, due to the interactions with the wake, effective 
forces on bubbles can even be non-local in time (history 
forces) [H, [H, \^ . Given these complications, for con- 
ceptional simplicity we decided to simply use the Auton 
lift coefficient Cl = 1/2 [4lj, realizing that at best qual- 
itative agreement between our simulations and possible 
experiments can be expected. Finally, periodic boundary 
conditions are employed. The aim of the present paper 
can thus only be to identify mechanisms; no quantitative 
predictions or comparisons with experiments are possi- 
ble. The numerical simulations for bubble-induced con- 
vection [2^ or for Taylor-Couette flow with microbubbles 
inducing drag reduction j42j were done in the same spirit 
with related numerical schemes, and indeed the relevant 
physical mechanisms could be identifled. 

Our numerical simulations are based on the code ex- 
tensively described in but now as stated above with 
forcing only by the bubbles, to mimick the pseudotur- 
bulent flow. For completeness we briefly repeat the dy- 
namical equations and the central assumption in section 
im Sections IIIII and IIVI describe the time evolution of 
global and spectral observables, respectively. In section 
fVl we propose a qualitative physical explanation for the 
detected fluid energy time evolution. Section IVTl contains 
conclusions. 



II. BUBBLES AND FLUID EQUATIONS 

A. Bubble motion equation 

The motion of a small, undeformable, bubble embed- 
ded in a velocity field u(x , t) can be modelled by the 
equation (see e.g. [J, [2^[4^): 

-^7=3— + -[u(y(i),i)-v(i)]-2g (2) 



~[v{t)^u{y{t),t)]xu;{y{t),t). 

The various symbols denote: y{t) the bubble position, 
V the bubble velocity, a; = V x u the fluid vorticity, 
g the gravity (acting in negative z-direction), and T(, the 
relaxation time, i.e. the time needed to adjust the bubble 
velocity to that of the fluid. The latter is related with the 
terminal velocity vt in still fluid, vt = '^grt- In the case 
of small bubble Reynolds number Re = 2|u — v|a/j^ < 1, 
with a the bubble radius and v the kinematic viscosity, 
for the bubble response time it holds Tb = a^/Gv 
IZsj . For larger bubbles the prefactor in this relation is 
larger than 1/6, which however would only quantitatively 
affect our results. The material derivative {D/Dt) of the 
fluid velocity is evaluated at the bubble position. Eq. 
^ embodies the effect of fluid acceleration plus added 
mass, drag force, buoyancy, and lift force, with the lift 
coefficient set to Cl = 1/2 for simplicity, as discussed 
in the introduction. We refer to refs. [ll, [13, [IB] for a 
detailed description of the terms in cq. 



B. Simulation of the flow 

The computational domain is a cube of side Lq = 2tt, 
consisting of — 128^ mesh points and subjected to 
periodic boundary conditions. The simulation is started 
at f = with the flow at rest and with Nt bubbles with 
Re ~ 0(1) placed at random locations. Bubbles rise 
because of gravity and transfer momentum to the fluid. 
We track their trajectories and we treat each bubble as a 
point-source of momentum. Then, the total action on the 
flow results by summin g th e (S-forcing that the bubbles 
apply at their positions [ITI. [26[. [43| : 

n / A \ 

(3) 

Here g is the gravity, directed along the negative z axis. 
The induced flow velocity u.{x, t) evolves according to the 
incompressible Navier-Stokes equation: 

u • Vu = -Vp + ^Au + ffc (4) 

which is solved by direct numerical simulation. 

We use the pseudo-spectral code described in [l3| 
where also the other details on the numerical simulations 
can be found. The point force approximation is validated 
by performing the same tests as in JJi]. We stress that, 
in contrast to that earlier work of ours, here eq. (|4|) does 
not contain any forcing on the large scales: the flow is 
sustained solely by the bubble forcing term f^. 

We analyse different cases. A list of the flow and bub- 
bles parameters is shown in Table [H In Table [IT] the 
values of the numerical parameters and their physical 
equivalents are presented for some of the simulations per- 
formed. 
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a 




Re 


Tb 


a (fim) 


Ep ■ 10^ 


E ■ 10 


(a) 


1.6% 


144000 


3.10 


8.4 • 10"^ 


78 


1.34 


1.96 


(b) 


0.8% 


72000 


3.10 


8.4 • 10"^ 


78 


0.67 


1.38 


(c) 


3.2% 


288000 


3.10 


8.4 ■ 10"^ 


78 


2.67 


2.25 


(d) 


1.6% 


50848 


4.37 


16.7 ■ 10^^ 


87 


1.34 


1.57 



TABLE I: Simulation parameters for all cases analysed: void 
fraction a, total bubbles number Ni,, bubble Reynolds num- 
ber Re = 2vTa./v, bubble response time Tb, equivalent radius 
a in physical units, potential flow estimate of the total energy 
induced in the flow Ep = av^/A, and time asymptotic esti- 
mate from our calculations of the total energy E. This value 
has been determined from time-averaging the statistically sta- 
tionary state. In the numerical simulations the kinematic vis- 
cosity u = 0.007, a, the rise speed in still fluid vt = 0.578, 
and Tb are fixed, the other quantities result consequently. In 
particular, the intensity of the gravity results from g = vr/^Tt 
and therefore is different in the simulation (d) as compared 
to the other ones. 





dimensionless 
parameter 


physical 
equivalent 


V 


0.007 


10"^cmVs 


g 


34.55 


981cm/ s"^ 


a 


0.019 


78/im 


Vt 


0.58 


2cm /s 


n 


8.4 • 10"=^ 


1ms 



TABLE IL Simulation parameters for cases (a) to (c) and cor- 
responding physical equivalents. The bubble Reynolds num- 
ber is Re — 2vTa/v = 3.10. 



III. EVOLUTION IN TIME FOR GLOBAL 
QUANTITIES 

We describe the energy time evolution of the pseudo- 
turbulent field generated by the rise of microbubbles. As 
already stated, at time t = 0, the bubbles are randomly 
placed in the flow, that is originally at rest. Their rise 
displaces liquid and thus generates velocity fluctuations. 
If these fluctuations are not rapidly dissipated by vis- 
cosity, they can be transmitted to larger scales. As a 
consequence, large scale motions are produced and the 
flow may become turbulent. 

We investigate this issue by measuring average flow 
quantities as well as by studying the spectral energy dis- 
tribution. We focus on case (a) of Table [D 

In Fig. [TJi we plot the total fluid energy, E{t) — 
{ux{tY + Uy{tY + Uz{t)'^)/2, as a function of time. In 
the beginning E{t) undergoes a steep rise, afterwards it 
slowly decreases, until it begins to oscillate and a statis- 
tically stationary state is reached. The behavior quali- 
tatively resembles that one of the front-tracking simula- 
tions of ref. 0, see figure 1 of that paper. The kinetic 
energy is mainly generated by the momentum transfer in 



the direction of gravity, as we confirm by measuring the 
three components of the fiuid velocity fiuctuations (uf), 
with i = x,y, z, which are much larger in the z direc- 
tion than in the horizontal ones, and the imidimensional 
Taylor-Reynolds number, defined as: 

Re\ = J ^^MI (no sum over i) . 

The behavior of as function of time is presented in 
Fig. [21 It is evident from the plot that the flow displays 
strong asymmetry. 

We now compare these results with potential flow the- 
ory, see e.g. the work of van Wijngaarden [1, IH, |4^. 
We do not expect agreement, given that we deal with 
point- par ticles and that the potential flow results of refs. 
[li IHj only hold in the high bubble Reynolds number 
case, and here we have Re between 3 and 5. Nonethe- 
less, we consider this comparison to be instructive and 
surprisingly we find the saturated kinetic energy to be of 
order aw|./4, in agreement with the potential flow theory 
results for high i?e bubbles. However, the redistribution 
of the energy along the three directions deviates from 
what is predicted by potential flow analysis, according to 
which we should have [isj : 

The potential flow value is {uj,) / {u^) ~ 4/3, whereas in 
our simulation this ratio is about 15. In the opposite 
limit. Re — *■ 0, under Stokes flow condition, the fluid 
equations are linear. In this regime, symmetry consider- 
ations require that rising bubbles cannot force the flow in 
directions perpendicular to their motion (provided that 
there are no walls but periodic boundary conditions as 
in our case), a result that is also intuitive for rising point 
particles in fluid at rest. As a consequence the ratio 
(uD/iUx) — *■ oo. Our result, for small but flnite Reynolds 
number, lies in between the two limits discussed. The 
same holds for the case of sedimenting particles in water 
(with larger particle Reynolds numbers between 38 and 
545), where this ratio is between 4 and 25 [30,] ■ 

The total energy induced by high Re bubbles, for which 
inertia effects are dominant, can be easily estimated by 
the following argument: at low void fractions the flow en- 
ergy is the sum of the energy induced by individual bub- 
bles, i.e., E ~ iVb(l/2)mbW^, where the effective mass 
of a bubble is mf, = p/(27ra'^/3), owing to the added 
mass factor 1/2. Thus E ~ av'^/A. On the other hand, 
when Re ~ 1, to estimate the total energy induced by 
sedimenting particles or rising bubbles is far more com- 
plicated. Indeed, for Re ^ 1, the flow induced by one 
particle decreases, with the distance r from the particle, 
as 1/r, thus leading to a total flow energy that diverges 
with the system size. Different screening mechanisms 
have been invoked in the past in order to account for 
this problem (see e.g. |5l|, [H, HI, IM III ) , but this 
issue is beyond the scope of the present paper. 
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FIG. 1: (a) Total fluid energy as function of time. The 
straight dashed line indicates the potential flow result: E — 
cra|n/4. (b) Total fluid energy in simulation with lift (solid 
line) and without (dashed line) as functions of time. 



In Table U we report the total energy estimated in our 
simulations, correspondent to different void fractions and 
bubble dimensions. In all cases the energy induced is of 
order av^/A. We note, by comparing cases (a) and (c?), 
that, when increasing the bubble dimension while fixing 
the void fraction, thus reducing the total bubble-fluid 
interface, less energy is generated in the flow. 



3.5 
3 

2.5 
2 

1.5 
1 

0.5 




3 10' 6 10' 9 10' 

t/x. 



FIG. 2: Behavior of the uni-dimensional Taylor-Reynolds 
number as function of time, in the a;-(pluses), {/-(circles), and 
2- (diamonds) directions. 
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IV. EVOLUTION IN TIME FOR SPECTRA 



After transforming to wavenumber space we consider 
the time development of the energy spectrum: 

E{k,t)^^ X! <(k,i)uj(k,t) i = x,y,z. (5) 

/£<|k|</£ + (//c 

Here Ui(k,t) is the ith component of the fluid velocity 
in fc-space and repeated indices are considered summed. 
E{k, t) is the total energy contained in a spherical shell 
of radius k and width dk. 

In Fig. [Ill the energy spectrum, averaged on four sub- 
sequent time intervals, is presented. The intervals corre- 
spond to < t/u < 6 • 10^, 12 • 10^ < t/n < 18 • 10^, 
24-10^ < t/n < 30-102, and lOO-lO^ < t/n < lOG-lO^ in 
Fig. [TJ The figure depicts that the energy is originally in- 
troduced at high wavenumbers and gradually transported 
to larger scales (solid line). However, after some time the 
spectrum flattens and a nearly constant energy is mea- 
sured at all scales. Thus, in the first stage of bubble- 
fluid coupling, an inverse energy cascade, from the small 



FIG. 3: (a) Energy spectra for the simulation that includes lift 
forces obtained by averaging on four different time intervals: 
< t/rt < 6-10^ (solid line), 12-10^ < t/n < 18-10^ (dashed 
line), 24 ■ 10^ < t/rb < 30 ■ 10^ (dotted line), and 100 • 10^ < 
t/Tb < 106 • 10^ (dot-dashed line). The straight line indicates 
the behavior in homogeneous and isotropic turbulence, (b) 
Energy spectra for the simulation without lift forces obtained 
by averaging on four different time intervals. For the various 
symbols look at the caption of Fig. [3ji. 



to the large scales, builds up large scale eddies. The 
corresponding slope of the spectrum - indicating the lo- 
cal transfer of energy - would be —5/3, however, as the 
Taylor- Reynolds numbers are small, such a scaling regime 
cannot really develop. Later on, the inverse energy cas- 
cade cannot be sustained and it flnally disappears in the 
flnal state. We investigate whether this last state is sta- 
tistically stationary by having a closer look at the energy 
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transfer equation in /c-space: 

^E{k, t) = T(k, t) - D{k, t) + Ffc(fc, t) (6) 
at 

Here, the various terms indicate, respectively: the energy 
transfer to wavenumber fc, T{k,t), 

T{k,t)^ J2 T{k,t), 

k<\'k\<k+dk 

where 

T(k,t) =Im ^/cju;'(k,i)^Uj(k-k',t)M,(k',t)^ , 

(7) 

the viscous dissipation, D{k,t) — 2i'k'^E{k,t), and the 
bubble forcing contribution, Fi,{k,t), 

Fb{k,t)= PbiKt), 

k<\k\<k+dk 

where 

Ffc(k,^) = 7^e(u*(k,^)•fb(k,^)), (8) 

with fb(k, t) the Fourier transform of the coupling term 
h{x, t), defined in eq. ([3|). Here Xm and TZe indicate, re- 
spectively, the imaginary and real part of the expression 
between the brackets. 

The spectra of the bubble forcing and of the dissipa- 
tion are shown in Fig. The strongest forcing is con- 
centrated on the small scales, as we expect, owing to 
the dimension of the energy sources. However, the en- 
ergy that is initially transfered to the large scales via the 
nonlinear interactions has to return to the small scales 
in order to be dissipated by viscosity. In fact, there is 
no energy sink on the large scales that could take the 
energy out of the system. Therefore the condition for 
establishing a stationary state is that the time average 
energy transfer has to be zero on all wavenumbers, and 
dissipation has to equal bubble forcing, i.e., T{k) = and 
D{k) = Fi,{k), where the time dependence has dropped 
out after the average on time. 

As we show in Fig. [4]2, apart from the large scales 
where the average still has not converged, this require- 
ment is satisfied by our simulation. In real flow, with 
walls instead of periodic boundary conditions, energy dis- 
sipation in the developing boundary layers will of course 
eventually play a crucial role in the energy balance. 

The time evolution of the energy spectrum can be com- 
pared to the one presented by 56], where the authors 
study a similar system, namely fluid motion generated by 
rising bubbles, by applying a different technique for the 
implementation of two-way coupling. The results agree 
qualitatively, i.e., the initial induction of structures at 
large scale is followed by a state in which the slope of the 
energy spectrum is reduced. As the authors state them- 
selves, the reason is connected with the temporal evolu- 
tion of the bubble distribution. Indeed, bubble clusters. 



that are assembled in the beginning and are able to force 
the liquid efficiently, are not stable and bubbles tend to 
distribute uniformly in the flow. Moreover, the struc- 
tures induced in the flow itself are far too weak to trap 
the bubbles. Thus, the phenomenon of vortex trapping 
of bubbles does not occur and therefore high local bub- 
ble concentration as in bubbly turbulent flows (see e.g. 
[l6, 57, 58]) are not created here. 

Note that the front-tracking simulations of ref. for 
27 bubbles with a bubble Reynolds number of about 25 
gave a quite different spectrum, namely a slope of about 
—3.6 in the large wavenumber regime. Indeed, a slope 
—3 has theoretically been attributed [9[ to this regime, 
in which the energy deposited by the wake is directly dis- 
sipated and which obviously cannot be modelled with the 
point-particle approach. All this wake energy dissipation 
is missing in our approach. 

We carry on the comparison with the results of [1^ , by 
looking at the high wavenumbers behavior of the spec- 
trum. In that paper a steeper slope with respect to ho- 
mogeneous and isotropic turbulence is observed. The 
critical wavenumber kc above which it shows up is es- 
timated by the average distance between the bubbles, 

1 /3 

that is Lc ^ 2n/N^' . In our simulation, for case (a) 
of Table |T1 we have: Lc ^ 0.12 (about 1/50 of the box 
with Lq = 27r), thus kc = 2tt /Lc ^ 52 and for case {d): 
Lc ^ 0.17, thus kc ~ 37. As we show in Fig. [51 a tran- 
sition in the slope of the energy spectrum occurs at high 
wavenumbers. However, neither the critical wavenumber 
nor the slope of the spectra can be clearly defined. 

We find agreement with the results of [5^ on the 
strongly anisotropic energy distribution along the three 
velocity components. Indeed, also in that work, about 
the 90% of the flow energy is contained in the vertical 
component (z) of the fluid velocity. 

We again also compare our results with the case for 
sedimenting particles in originally still fluid [13, HH : Also 
for that case the energy spectrum is very broad. The 
spectral slope of the frequency spectrum was found to 
be around —1.1 for remarkably several decades, but a 
comparison of this value with the slopes in the wavenum- 
ber spectra obtained in this paper is difhcult as Taylor's 
frozen flow hypothesis will most likely not hold for this 
weak turbulence. 



V. PHYSICAL EXPLANATION OF THE 
RESULTS 

The occurrence of the inverse cascade phenomenon in 
three dimensional turbulence has been related to the 
presence of strong anisotropics at small scales (see e.g. 
[59, 60]). These anisotropics can be produced by bub- 
ble clusters elongated in the gravity direction. Within 
them, the energy production term due to the bubbles 
can be far more intense in the vertical direction than in 
the horizontal ones, owing to the high values reached by 
the (u • g) contribution. The stability of these structures 
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FIG. 4: (a) Time average of the contribution of the bubble 
forcing to the energy spectrum, as defined in eq. ((8)1 (soHd 
Une) and of the viscous energy dissipation D{k) = 2uk^E{k) 
(dotted Une), in the simulation with lift force. (6) Time av- 
erage of the contribution of the bubble forcing to the energy 
spectrum (solid line) and of the viscous energy dissipation 
D{k) — 2vk^ E{k) (dotted line), in the simulation without lift 
force. 



is opposed by horizontal forces that laterally spread the 
bubbles. When considering the bubble motion equation, 
it appears that the lift is the most relevant of such forces. 
Therefore we further investigate the system by compar- 
ing the outcome of simulations with lift and without lift 
force. 

Surprisingly, without lift the results are very different. 
In Fig. [H) the total energy in the two simulations, one 
including the lift (solid line) and the other excluding it 
(dashed line), are compared. In both cases the bubbles 
parameters correspond to run (a) in Table HI The energy 
induced in the second simulation is up to 30 times larger 
than in the first. Also Murai et al. [61| measured a higher 
turbulence intensity in numerical simulations without lift 
force than in simulations with lift, though the difference 
detected is quantitatively much smaller than in our case. 

The behavior in spectral space is remarkably different, 
too. In Fig.[3f) the time evolution of the energy spectrum 
in the latter simulation is presented. The spectrum is 
averaged on four subsequent time intervals, which are the 
same as for FiglSli. The process of inverse energy cascade 
is now strongly enhanced. In fact, the spectral intensity 
at small wavenumbers increases in time, whereas it is 
constant at large ones. 

Moreover, it is remarkable that a local slope close to 
—5/3 settles nearly at once at high wavenumbers and 
is stable during the whole process - though the scaling 




FIG. 5: Energy spectra, in the statistically stationary state, 
for case (a) (dashed line) and (d) (solid line) of Table IH 



regime again is not very pronounced due to the small 
Taylor-Reynolds numbers. Therefore the small scale forc- 
ing is strong enough to generate a flow that presents the 
same characteristics as real three-dimensional turbulence 
- or also two-dimensional turbulence in the inverse cas- 
cade regime. On the other hand, this simulation is not 
statistically stationary. In Fig. we show that there is 
a difference on a wide range of scales between the bub- 
bles forcing term Fh{k) and the fluid viscous dissipation 
D{k). Thus the condition of stationarity is not fulfilled 
and the large flow scales are still fed with energy from 
the small ones. 

We again stress that the model for the equation 
of motion without lift force does not give a complete 
representation of the surface forces acting on bubbles 
or particles of Re ~ 0{1). Indeed, previous works (see 
e.g. refs. [13, Hlj and ref. [i^l for the case of uniform 
shear) have pointed out the relevance of the lift in this 
regime. Of course our results can only expected to be 
qualitative and not quantitative, as the near-field inter- 
actions between the bubbles is not correctly described 
by the present point-particle approach and as the lift 
coefficient is set to be constant, Cl = 1/2. However, 
refined expressions for the lift force arc not likely to give 
qualitatively different behaviors. The main effect of the 
lift is to cause the bubble dispersion along the horizontal 
directions, thus strongly reducing the anisotropy in the 
flow caused by the forcing in the vertical direction. 
Indeed, by definition, the hft force drifts the bubbles in 
horizontal planes, in directions perpendicular to their 
average motion. 

We now qualitatively investigate the breaking effect of 
the lift on vertical bubble chains. Note again that for a 
quantitative analysis the near-field as e.g. obtained from 
the front-tracking simulations of refs. [y, 0, [H, is 
crucial. We base our analysis on the description of the 
two-bubble long range interactions proposed in ref. [50| . 
The main finding is that a bubble in the wake of another 
one experiences, because of the lift, a lateral force, lead- 
ing to a deficit of nearby bubbles in the gravity direction 
(see Fig.lH). 

We further explore this phenomenon by computing the 
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bubble 
trajectory 





FIG. 6: Sketch of the action of the hft force on a bubble rising 
in the wake of another one: the hft tends to expel the bubble 
from the wake. 



bubble density autocorrelation function, defined accord- 
ing to: 



P12(R) = 



(c^(x + RK(x)) 
(c'(x)c'(x)) 



(9) 



Here c'(x) is the fluctuation of the bubble concentration 
in X with respect to the average value a, and the brackets 
denote averages over all x. We consider the autocorre- 
lation in the horizontal, (x-y), plane and in the vertical, 
(z), direction separately. The analysis is carried out for 
the two simulations presented in Fig. [2). The results 
are plotted in Fig. [T] It is shown that, whereas in the 
simulation without lift forces the pair correlation goes 
monotonically to zero, in the one with lift the autocor- 
relation in the z-direction becomes negative and later on 
approaches zero from below. A negative autocorrelation 
at small distances R is detected also in the horizontal 
planes. The interpretation of the result is that the bub- 
bles approach is resisted by the lift, and this is occurring 
especially in the vertical direction, where a bubble rising 
in the wake of another experiences horizontal forces that 
expels it from the wake. 

Qualitatively, the organization of the bubbles in our 
simplified simulations is similar to what has been ob- 
served in above mentioned front tracking simulations by 
Tryggvason and coworkers [1, 0, [13, [H, [2§| or in the ex- 
periments of ref. [13]: Small bubbles with Re ~ 1 were 
dispersed in a nearly homogeneous manner, with an in- 
creasing tendency of horizontal alignment when the bub- 
ble Reynolds number approaches 10. Also in the recent 
experiments by Harteveld et al. (63j bubbly driven flows 
have been found to be rather homogeneous and no vortex 
trapping has been detected. The front-tracking simula- 
tions show that for even larger bubble Reynolds num- 
bers horizontal bubble clusters can emerge, which also 
have been seen in the experiments of [25j . Even in the 
two-dimensional simulations described in ref. (6^ a pref- 
erential tendency of bubbles to stay rather side-by-side 
(along the horizontal direction) than in "tandem config- 



FIG. 7: Autocorrelation function pi2(R) as a function of the 
distance |R|, in the horizontal x-y directions (open symbols) 
and in the vertical z direction (filled symbols). The results 
indicate simulations without hft force (diamonds) and with 
hft force (squares). 



uration" (along the vertical) was reported. 



VI. CONCLUSIONS 

The behavior of a flow driven exclusively by rising bub- 
bles has been investigated by direct numerical simulation 
for the Navier Stokes equations and Lagrangian tracking 
for the bubble trajectories, where the bubbles have been 
treated as point particles on which effective forces act. 
The evolution of global quantities, like the total flow en- 
ergy, as well as of spectral quantities has been followed in 
time. The results show that the bubble motion initially 
generates large scale structures by local in scale energy 
transfer, though the corresponding —5/3 scaling regime 
in the spectrum is not very pronounced due to the small 
Taylor-Reynolds numbers. Later on, however, the bub- 
bles distribution tends to be more disperse, the energy 
spectrum becomes flat and energy input equal viscous 
dissipation at all scales. Therefore, the statistically sta- 
tionary state of this pseudo-turbulent velocity field does 
not possess the characteristics of real turbulence. 

We give qualitative evidence that the physics that de- 
termines it are the bubble-bubble indirect interactions 
that occur via the carrier flow. Indeed, a bubble in the 
wake of another one, experiences, because of the lift, a 
horizontal force that prevents the assembling and the sta- 
bility of vertical clusters, similarity as has been observed 
in the front-tracking simulations of bubbles of compara- 
ble size [13, [13 or of larger bubbles [1, [3| , where the bub- 
bles also distribute homogeneously and horizontal pairs 
of bubbles are favored. In our case, as a consequence the 
total forcing induced in the flow is not strong enough to 
sustain high energy levels and an inverse energy cascade 
from small to large scales. 

The results presented in this work apply to a flow with 
periodic boundary conditions. In a real experiment such 
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as in Taylor-Couette flow, the existence of boundaries can 
lead to the generation of large scale vortex structures, 
that, in turn, affect the bubble motion, as seen in the 
experiments of Murai et al. 65*1 and the corresponding 
numerical simulations of Sugiyama et al. [42|. In those 
flows the energy dissipation in the developing boundary 
layers will play a crucial role in the stationary energy 
balance. 

We understand this work as complementary to the 
front-tracking simulations: Bunner and Tryggvason [6] 
end their numerical study on the dynamics of homoge- 
neous bubble flows with the question: "What happens 
when the number of bubbles increases beyond 216?" (the 
number of bubbles they could numerically treat in the pa- 



per). Here we treat up to 288000 bubbles, though in an 
approximate and simplifled way. Nonetheless, we see very 
similar phenomena as observed in the front-tracking sim- 
ulations, and can even identify the lift force as origin of 
the homogeneous bubble distribution, by artifically turn- 
ing it off. A full verification can of course only come from 
simulations which both resolve the individual bubbles in- 
cluding their wakes and deal with hundred-thousands of 
individual bubbles. Such simulations will however unfor- 
tunately not be numerically feasible for many years to 
come. 
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